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(54) Title: SPECTRAL SUBTRACTION NOISE SUPPRESSION METHOD 
(57) Abstract 

A spectral subtraction noise suppression method in a frame 
based digital communication system is described. Each frame includes 
a predetermined number N of audio samples, thereby giving each 
frame N degrees of freedom. The method is performed by a spectral 
subtraction (150) function #(w) which is based on an estimate (140) 
$v(w) of the power spectral density of background noise of non-speech 
frames and an estimate (130) *»(w) of the power spectral density 
of speech frames. Each speech frame is approximated (120) by a 
parametric model that reduces the number of degrees of freedom to 
less than N. The estimate $*(*) of the power spectral density of each 
speech frame is estimated (130) from the approximative parametric 
model. 
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SPECTRAL SUBTRACTION NOISE SUPPRESSION METHOD 

TECHNICAL FIELD 
The present invention relates to noise suppresion in digita! frame based communication 
systems, and ,n particular to a spectral subtraction noise suppression method ,n such 

systems. 

BACKGROUND OF THE INVENTION 
A common problem in speech signal processing is the enhancement of a speech signal 
from its noisy measurement. One approach for speech enhancement based on angle 
channel (microphone) measurements is filtering in the frequency domain applying spectra! 
subtraction techniques, W . |2|. Under the assumption that the background noise is long- 
time stationary (in compare with the speech) a mode, of the background noise is usually 
esfmated duringtime interval with non-speech activity. Then, during data frames with 
speech activity, this estimated noise model is used together with an estimated model of 
' the noisy speech in order to enhance the speech. For the spectral subtraction techmques 
these models are traditional* gwen in terms of the Power Spectral Density (PSD), that 
is estimated using classical FFT methods. 

None of the abovementioned techniques give in their basic form an output signal with 
satisfactory audible quality in mobile" telephony applications, that is 

1. non distorted speech output 

2. sufficient, reduction of the noise level 

3. remaining noise without annoying artifacts 

,„ particular, the spectra, subtraction methods are known to vioiate > when 2 is fulfilled 
or v,o,.te 2 when . is fulfilled. In addition, in most cases 3 is more or less violated su.ce 
the methods introduce, so called, musical noise. 

The above drawbacks with the spectral subtraction methods have been known and, 
in the literature, several ad hoc modifications of the basic algorithms have appeared 
for particular speech-in-noise scenarios. However, the problem how to design a spectra, 
subtraction method that for general scenarios fulfills 1-3 has remained unsoived. . 
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In order to highlight the difficulties with speech enhancement from noisy data, note 
that the spectral subtraction methods are based on filtering using estimated models of the 
incoming data. If those estimated models are close to the underlying "true" models this 
is a well working approach. However, due to the short time stationary of the speech (10- 
40 ms) as well as the physical reality surrounding a mobile telephony application (8000Hz 
sampling frequency, 0.5-2.0 s stationarity of the noise, etc.)' the estimated models are 
likely to significantly differ from the underlying reality and, thus, result in a filtered 
output with low audible quality. 

EP, Al, 0 588 526 describes a method in which spectral analysis is performed either 
with Fast Fourier Transformation (FFT) or Linear Predictive Coding (LPG). 

SUMMARY OF THE INVENTION 

An object of the present invention is to provide a spectral subtraction noise suppresion 
method that gives a better noise reduction without sacrificing audible quality. 
This object is solved by the characterizing features of claim 1. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The invention, together with further objects and advantages thereof, may best be 
understood by making reference to the following description taken together with the 
accompanying drawings, in which: - 

FIGURE 1 is a block diagram of a spectral subtraction noise suppression system 
suitable for performing the method of the present invention; 

FIGURE 2 is a state diagram of a Voice Activity Detector (VAD) that may be used 
in the system of Fig. 1; 

FIGURE 3 is a diagram of two different Power Spectrum Density estimates of a speech 
frame; 

FIGURE 4 is a time diagram of a sampled audio signal containing speech and back- 
ground noise; 

FIGURE 5 is a time diagram of the signal in Fig. 3 after spectral noise subtraction 
in accordance with the prior art; 

FIGURE 6 is a time diagram of the signal in Fig. 3 after spectral noise subtraction 
in accordance with the present invention; and 

FIGURE 7 is a flow chart illustrating the method of the present invention. 
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DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

THE SPECTRAL SUBTRACTION TECHNIQUE 
Consider a frame of speech degraded by additive noise 

il» -.<*) + »(*) *-l W (1) 

where *(*), .(*) and denote, respectively, the noisy measurement of the speech, the 
speech and the additive noise, and N denotes the number of samples in a frame. 

The speech is assumed stationary over the frame, while the noise is assumed long- 
time stationary, that is stationary over several frames. The number of frames where » W 
is stationary is denoted by x » 1. Father, it is assumed that the speech activity ,s 
sufficientlylow.sothatamodelofthenoiseca.beaccuratelyestimatedduringnon-speech 

activity. 

i j f*p<;n«^ of respectively, the measurement, tne 

Denote the power spectral densities (PSDs) oi, respecu j, 

speech and the noise by ♦,(«), *.(«) and •„(«-). where 

»,(w) = *.H + *» (2) 

Knowin B ♦,(«) and ♦.(„). the quantities •.(«) and . (*) can be estimated using standard 
spectral subtraction methods, cf [2], shortly reviewed below 
Let s{k) denote an estimate of s{k). Then, 

(3) 

where W denotes some linear transform, for example the Discrete Fourier Transform 
(DPT) and where *i» is a real-valued even function in „ 6 (0,2.) and such that. 
0 < HW) < I- The function if (*) depends on ♦.<-) and ♦.(«). Since if < w ) is real- 
valued, the phase of 5(.) = if H *M equals the phase of the degraded speech. The use 
of real-valued if ( W ) is motivated by the human ears unsensitivity for phase distortion. 

In general, ♦.(„) and *M are unknown and have to be replaced in *(,) by est,- 
mated quantities 4.H and *.(«). Due to the non-stationarity of the speech, » 
estimated from a single frame of data, while ♦.(„) is estimated using data in r speech free 
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frames. For simplicity, it is assumed that a Voice Activity Detector (VAD) is available in 
order to distinguish between frames containing noisy speech and frames containing noise 
only. It is assumed that $ v (lj) is estimated during non-speech activity by averaging over 
several frames, for example, using 

= P* v (u) e ~ l + (1 - p)* v (w) (4) 

In (4), $ v (w)' is the (running) averaged PSD estimate based on data up to and including 
frame number I and $ v (u») is the estimate based on the current frame. The scalar p € (0, 1) 
is tuned in relation to the assumed stationary of v(Ar). An average over r frames roughly 
corresponds to p implicitly given by 

A suitable PSD estimate (assuming no apriori assumptions on the spectral shape of the 



*.M = jfV( u )V(u) ( 6 ) 

where "*" denotes the complex conjugate and where V(u) = P(v(k)). With F{ ) = 
FFT(-) (Fast Fourier Transformation), $„(w) is the Periodigram and 4>„(w) in (4) is the 
averaged Periodigram, both leading to asymptotically (N » 1) unbiased PSD estimates 
with approximative variances 



Var(* u H) » 
Var(*.(w)) * 

T 

A similar expression to (7) holds true for 4> x (u;) during speech activity (replacing $>l(u>) 
in (7) with *2(u;))- 

A spectral subtraction noise suppression system suitable for performing the method 
of the present invention is illustrated in block form in Fig. 1. From a microphone 10 
the audio signal x(t) is forwarded to an A/D converter 12. A/D converter 12 forwards 
digitized audio samples in frame form {x(k)} to a transform block 14, for example a 
FFT (Fast Fourier Transform) block, which transforms each frame into a corresponding 
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freq «encv transformed frame (X<*)>. The transformed frame is filtered by A( w ) in block 
16 This step performs the actual spectral subtraction. The resulting signal { S(u)} is 
transformed back to the time domain by an inverse transform block 18. The result is a 
frame <S(*.)> in which the noise has been suppressed. This frame may be forwarded to 
an echo canceler 20 and thereafter to a speech encoder 22. The speech encoded signal is 
then forwarded to a channel encoder and modulator for transmission (these elements are 
not shown). 

The actual form of H{u) in block 16 depends on the estimates •,(«), •.(«), which 
are formed in PSD estimator 24, and the analytical expression of these estimates that 
is used. Examples of different expressions are given in Table 2 of the next section. The 
major part of the following description will concentrate on different methods of fanning 
estimates * x (w), *» from the input frame {*(*)}• 

PSD estimator 24 is controlled by a Voice Activity Detector (VAD) 26, which u SeS 
input frame {x(fc» to determine whether the frame contains speech (S) or background 
noise (B). A suitable VAD is described in (5), [6]. The VAD may be implemented as 
a state machine having the 4 states illustrated in Fig. 2. The resulting control signal 
S/B is forwarded to PSD estimator 24. When VAD 26 indicates speech (S), states 21 
and 22 PSD estimator 24 will form *.(«). On the other hand, when VAD 26 indicates 
non-speech activity (B). state 20, PSD estimator 24 will form *.(«). The latter estimate 
will be used to form H{«) during the next speech frame sequence (together w.th *,<w) 
of each of the frames of that sequence). 

Signal S/B is also forwarded to spectral subtraction block 16. In this way block 16 
mav Ipply different filters during speech and. non-speech frames. During speech frames 
HW) is the above mentioned expression of *,<»). *.(«). On the other hand, dunng 
non-speech frames *(«> may be a constant H (0 < H < D that reduces the background 
sound level to the same level as the background sound level that remains in speech frames 

t„ thi, wav the perceived noise level will be the same during both 
after noise suppression. In this way tne perceivcu 

speech and non-speech frames. 

Before the output signal J(t) in (3) is calcuiated, 6M »•* >» • emb ° d '- 
ment, be post filtered according to 

• /r,(w).nax(fl.l,W(«W) Vw (8) 
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Table 1: The postfiltering functions. 



STATE (st) 


H(u) 


COMMENT 


0 


1 (Vu/) 


s(k) = x(k) 


20 


0.316 (Vw) 


muting -lOdB 


21 


0.7 H(u) 


cautios filtering (-3dB) 


22 







where H{u) is calculated according to Table 1. The scalar 0.1 implies that the noise floor 
is -20dB. 

Furthermore, signal S/B is also forwarded to speech encoder 22. This enables different 
encoding of speech and background sounds. 

PSD ERROR ANALYSIS 

It is obvious that the stationarity assumptions imposed on s(k) and v(k) give rise to 
bounds on how accurate the estimate s(k) is in comparison with the noise free speech signal 
s(k). In this Section, an analysis technique for spectral subtraction methods is introduced. 
It is based on first order approximations of the PSD estimates $ x (u>) and, respectively, 
$ v (<jj) (see (11) below), in combination with approximative (zero order approximations) 
expressions for the accuracy of the introduced deviations. Explicitly, in the following an 
expression is derived for the frequency domain error of the estimated signal s(A:), due to 
the method used (the choice of transfer function H{u)) and due to the accuracy of the 
involved PSD estimators. Due to the human ears unsensitivity for phase distortion it is 
relevant to consider the PSD error, defined by 

= *«(") - *.M (9) 

where 

$ s (u>) = H 7 (u)*:(") (10) 

Note that $ 3 (v) by construction is an error term describing the difference (in the frequency 
domain) between the magnitude of the filtered noisy measurement and the magnitude of 
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the speech. Therefore, ♦.(«) can take both positive and negative values and is not the 
PSD of any time domain signal. In (10), denotes an estimate of H{u) based on *,(«) 
and In this Section, the analysis is restricted to the case of Power Subtraction 

(PS), [2]. Other choices of (w) can be analyzed in a similar way (see APPENDIX A-C). 
In addition novel choices of H{*>) are introduced and analyzed (see APPENDIX D-G). A 
summary of different suitable choices of H{w) is given in Table 2. 

Table 2: Examples of different spectral subtraction methods: Power Sub- 
traction (PS) (standard PS, HpsI*>) *» 6 = *)' Magnitude Sub- 
traction (MS), spectral subtraction methods based on Wiener Fil- 
tering (WF) and Maximum Likelihood (ML) methodologies and 
Improved Power Subtraction (IPS) in accordance with a preferred 
embodiment of the present invention. 

H(«) 



H 6PS [u>) = - «*„(w)/*x(") 



Hwf{") = 



By definition, //(c) belongs to the interval 0 < H(u>) < 1, which not. necesarilly 
holds true for the corresponding estimated quantities in Table 2 and, therfore, in practice 
half-wave or full-wave rectification, (1], is used. 

In order to perform the analysis, assume that the frame length N is sufficiently large 
(N » 1) so that ♦,(«) and *.(«) are approximately unbiased. Introduce the first order 
deviations 

* x (u/) = + A r (w) 

(11) 
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**(") = *«(w) + A v (w) 

where A x (u/) and A v (u) are zero-mean stochastic variables such that 
£[A x (u/)/#,(<j)p « i and E[A v (u)/$ v (uj)) 2 « 1. Here and in the sequel, the notation 
E[] denotes statistical expectation. Further, if the correlation time of the noise is short 
compared to the frame length, £[(*„(*)'- - # w (u;))] % 0 for ^ fc, where 

is the estimate based on the data in the Mh frame. This implies that A r (u/) 
and A v (u>) are approximately independent. Otherwise, if the noise is strongly correlated, 
assume that $ v {u) has a limited (<£ N) number of (strong) peaks located at frequencies 

wi w„. Then, £[(* w (u/)'- (*„((*;)* - * tf (u/))] % 0 holds for w ^ w y j = 1, . . . , n 

and f 7^ A: and the analysis still holds true for u) ^ u s j = 1, . . . , n. 

Equation (11) implies that .asymptotical (AT > l) unbiased PSD estimators such as 
the Periodogram or the averaged Periodogram are used. However, using asymptotically 
biased PSD estimators, such as the Blackman- Turkey PSD estimator, a similar analysis 
holds true replacing (11) with 

*x(w) = *,(u/) + A x (w) + B x (u) 

and 



= *„(u;) + A B (u/) + B„(u;) 



where, respectively, B x (u>) and 5 v (u;) are deterministic terms describing the asymptotic 
bias in the PSD estimators. 

Further, equation (11) implies that $ 3 [u) in (9) is (in the first order approximation) 
a linear function in A x ((j) and A v (w). In the following, the performance of the different 
methods in terms of the bias error (£[$ 3 (u/)]) and the error variance (Var(<t> s (u;))) are 
considered. A complete derivation will be given for H ps (uj) in the next section. Similar 
derivations for the other spectral subtraction methods of Table 1 are given in APPENDIX 
A-G. 

ANALYSIS OF H PS (u>) (H 6 ps(") for 6=1) 

Inserting (10) and Hps{v) from Table 2 into (9). using the Taylor series expansion 
(Hi)" 1 ^ l-i and neglecting higher than first order deviations, a straightforward 
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calculation gives 

where is used to denote an approximate equality in which only the dominant terms 
are retained. The quantities M") and A u (u,) are zero-mean stochastic variables. Thus, 

and 

Var^M) * £fj$V.r(*,M) + Var (*^» ( 14 > 
In order to continue we use the general result that, for an asymptotically unbiased spectral 
estimator *(u>), cf (7) 

Var($(u»))=:7(w)* 2 (w) ( 15 ) 

for some (possibly frequency dependent) variable 7 (")- For example, the Periodogram 
corresponds to 7 ( W ) « 1 + (sin** /Nsinu)\ which for N » 1 reduces to 7 » 1. 
Combining (14) and (15) gives 

" Var(4.(u)) = 7 (l6) 



RESULTS FOR HmsM 

Similar calculations for HmsM give (details are given in APPENDIX A): 



and 



Var(*.(*)) a ll- V 1 + $j^) 



2 



1 + 12^1 7 *»(u>) 



RESULTS FOR Hwf(v) 

Calculations for HwfH give (details are given in 
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*v(") 



and 



Var(i. M ).4(l-|^|) 2 7 ^ (u;) 



RESULTS FOR H ML {u>) 

Calculations for H ml {uj) give (details are given in APPENDIX C): 



and 



Var($,( w )) - ^ ^ + ^ 



RESULTS FOR H lPS {u) 

Calculations for H lPS (u) give (#/ PS (u;) is derived in APPENDIX D and analyzed 
APPENDIX E): 



£[*.(«)! =:(<?(«) -!)*,(<</) 



and 



Var(* f (u;)) a <7 2 (u,) 



COMMON FEATURES 
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For the considered methods it is noted that the bias error only depends on the choice 
of *(*). while the error variance depends both on the choice of Jf( w ) and the variance of 
the PSD estimators used. For example, for the averaged Periodogram estimate of •„(«) 
one has, from (7), that lv « 1/r. On the other hand, using a single frame Periodogram 
for the estimation of ♦.(*), one has 7, « 1- Thus, for r » 1 the dominant term in 
T - 7, + 7v, appearing in the above vriance equations, is 7x and thus the main error 
source's the single frame PSD estimate based on the the noisy speech. 

Prom the above remarks, it follows that in order to improve the spectral subtraction 
techniques, it is desirable to decrease the value of 7 . (select an appropriate PSD estimator, 
that is an approximately unbiased estimator with as good performance as possible) and 
select a "good" spectral subtraction technique (select *<«)). A key idea of the present 
invention is that the value of lx can be reduced using physical modeling (reducing the 
number of degrees of freedom from N (the number of samples in a frame) to a value less 
than N) of the vocal tract. It is well known that s(k) can be accurately described by an 
autoregressive (AR) model (typically of order P « 10). This is the topic of the next two 
sections. 

In addition, the accuracy of *.(*) (and, implicitly, the accuracy of s(k)) depend, 
on the choice of H[u). New, preferred choices of H(«) are derived and analyzed a 
APPENDIX D-G. 

SPEECH AR MODELING 

In a preferred embodiment of the present invention .(*)!. modeled as an autoregressive 

(AR) process 

s(l) = _L_. (k) N .(17) 

where *(,-') is a monie (the leading coefficient equals one) n-th order P olynom,al in the 
backward shift operator fa'M*) = »<* - ] )' elc \ 

/!(,-') = 1 + «W"' +••■ + V?" P (18> 
and «,(*) is white zero-mean noise with variance „>. At a first g.a»ce, it may seem re- 
strictive to consider AR models only. However, the use of AR models for speech modeling 
is motivated both from physical modeling of the vocal tract and, which is more important 
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here, from physical limitations from the noisy speech on the accuracy of rhe estimated 
models. 

In speech signal processing, the frame length N may not be large enough to allow 
application of averaging techniques inside the frame in order to reduce the variance and, 
still, preserve the unbiasness of the PSD estimator. Thus, in order to decrease the effect 
of the first term in for example equation (12) physical modeling of the vocal tract has to 
be used. The AR structure (17) is imposed onto s(k). Explicitly, 

*' M = M(^)j5 + *»M .(») 
In addition, $ v (u) may be described with a parametric model 

where B(q' l ) : and C(q~ l ) are, respectively, g-th and r-th order polynomials, defined 
similarly to A(q~ l ) in (18). For simplicity a parametric noise model in (20) is used in 
the discussion below where the order of the parametric model is estimated. However, it 
is appreciated that other models of background noise are also possible. Combining (19) 
and (20), one can show that 

where i)(k) is zero mean white noise with variance a\ and where D{q~ x ) is given by the 
identity 

. ^ 2 |Z>(0| 2 = a 2 w \C(e^)\ 2 ■+ *l\B(e")\ 2 \A(e")\ 2 (22) 

SPEECH PARAMETER ESTIMATION 

Estimating the parameters in (17)-(18) is straightforward when no additional noise is 
present. Note that in the noise free case, the second term on the right hand side of (22) 
vanishes and, thus, (21) reduces to (17) after pole-zero cancellations. 

Here, a PSD estimator based on the autocorrelation method is sought. The motivation 
for this is fourfold. 



The autocorrelation method is well known. In particular, the estimated parameters 
are minimum phase, ensuring the stability of the resulting filter. 
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. Using the Levinson algorithm, the method is easily implemented and has a low 

computational complexity. 
. An optimal procedure includes a nonlinear optimization, explicitly requiring some 

initialization procedure. The autocorrelation method requires none. 

. . From a practical point of view, it is favorable if the same estimation procedure 
can be used for the degraded speech and, respectively, the clean speech when it 
is available. In other words, the estimation method should be independent of the 
actual scenario of operation, that is independent of the speech-to-noise ratio. 

It is well known that an ARMA model (such as (21)) can be modeled by an infinite 
order AR process. When a finite number of data are available for parameter estimation, 
the infinite order AR model has to be truncated. Here, the model used is 

*>-7(PJ«W (23) 

where F(g-') is of order p. An appropriate model order follows from the discussion below. 
The approximative model (23) is close to the speech in noise process if their PSDs are 
approximately equal, that is 

\D(e^)\ 2 „ 1 (24) 
|>l(e-)| 2 |C(e-)| 2 ~\|F(e-)| 2 

Based on the physical modeling of. the vocal tract, it is common to consider p = 
deg(A(<7- 1 )) = 10. From (24) it also follows that p = deg(F( 9 " 1 ) » deg(A(<T 1 )) + 
deg(C( 7 - 1 )) = P + r, where p + r roughly equals the number of peaks in ♦,(<*)■ On the 
other hand, modeling noisy narrow band processes using AR models requires f « N in 
order to ensure realible PSD estimates. Summarizing, 

p + r 

A suitable rule-of-thumb is given by p ~ JR. From the above discussion, one can expect 
that a parametric approach is fruitful when N » 100. One can also conclude from (22) 
that the flatter the noise spectra is the smaller values of N is allowed. Even if p is not 
large enough, the parametric approach is expected to give reasonable results. The reason 
for this is that the parametric approach gives, in terms of error variance, significantly 
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more accurate PSD estimates than a Periodogram based approach (in a typical example 
the ratio between the variances equals 1:8; see below), which significantly reduce artifacts 
as tonal noise in the output. 

The parametric PSD estimator is summarized as follows. Use the autocorrelation 
method and a high order AR model (model order p » p and p ~ t/n) in order to 
calculate the AR parameters {fi>...J p } and the noise variance a\ in (23). From the 
estimated AR model calculate (in N discrete points corresponding to the frequency bins 
of X(u) in (3)) 4> z (u>) according to 

*' M = i7& < 25 > 

Then one of the considered spectral subtraction techniques in Table 2 is used in order to 
enhance the speech $(k). 

Next a low order approximation for the variance of the parametric PSD estimator 
(similar to (7) for the nonparametric methods considered) and, thus, a Fourier series ex- 
pansion of s(k) is used under the assumption that the noise is white. Then the asymptotic 
(for both the number of data (N » 1) and the model order (p » 1)) variance of $ x (u;) 
is given by 

Var(*«H) |^(u,) (26) 

The above expression also holds true for a pure (high-order) AR process. From (26), it 
directly follows that 7x % 2p/h\ that, according to the aforementioned rule-of-thumb, 
approximately equals 7x - which should be compared with 7* as 1 that holds true 

for a Periodogram based PSD estimator. 

As an example, in a mobile telephony hands free environment, it is reasonable to 
assume that the noise is stationary for about 0.5 s (at 8000 Hz sampling rate and frame 
length N = 256) that gives r % 15 and, thus, j v ~ 1/15. Further, for p = y/N we have 
7, = 1/8. 

Fig. 3 illustrates the difference between a periodogram PSD estimate and a parametric 
PSD estimate in accordance with the present invention for a typical speech frame. In this 
example N=256 (256 samples) and an AR model with 10 parameters has been used. It is 
noted that the parametric PSD estimate $ x (u) is much smoother than the corresponding 
periodogram PSD estimate. 
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Fig. 4 illustrates 5 seconds of a sampled audio signal containing speech in a noisy 
background. Fig. 5 illustrates the signal of Fig. 4 after spectral subtraction based on a 
periodogram PSD estimate that gives priority to high audible quality. Fig. 6 illustrates 
the signal of Fig. 4 after spectral subtraction based on a parametric PSD estimate in 
accordance with the present invention. 

A comparison of Fig. 5 and Fig. 6 shows that a significant noise suppression (of the 
order of 10 dB) is obtained by the method in accordance with the present invention. (As 
was noted above in connection with the description of Fig. 1 the reduced noise levels 
are the same in both speech and non-speech frames.) Another difference, which is not 
apparent from Fig.' 6, is that the resulting speech signal is less distorted than the speech 
signal of Fig. 5. 

The theoretical results, in terms of bias and error variance of the PSD error, for all 
the considered methods are summarized in Table 3. 

It is possible to rank the different methods. One can, at least, distinguish two criteria 
for how to select an appropriate method. 

.First, for low instantaneous SNR, it is desirable that the method has low variance in 
order to avoid tonal artifacts in s(Jfc). This is not possible without an increased bias, and 
this bias term should, in order to suppress (and not amplify) the frequency regions with 
low instantaneous SNR, have a negative sign (thus, forcing 4>» in (9) towards zero). 
The candidates that fulfill this criterion are, respectively, MS, IPS and WF. 

Secondly, for high instantaneous SNR, a low rate of speech distortion is desirable. 
Further if the bias term is dominant, it should have a positive sign. ML, 6PS, PS, IPS 
and (possibly) WF fulfill the first statement. The bias term dominates in the MSE 
expression only for ML and WF, where the sign of the bias terms are positive for ML 
and, respectively, negative for WF. Thus, ML, 6PS, PS and IPS fulfill this criterion. 

ALGORITHMIC ASPECTS 

In this section preferred embodiments of the spectral subtraction method in accordance 
with the present invention are described with reference to Fig. 7. 



1. Input: x= {x(k)\k = 1 N}. 

2. Design variables 



IG 
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Table 3: Bias and variance expressions for Power Subtraction (PS) (stan- 
dard PS, Hpsfa) for 6 = 1), Magnitude subtraction (MS), Im- 
proved Power Subtraction (IPS) and spectral subtraction meth- 
ods based on Wiener Filtering (WF) and Maximum Likelihood 
(ML) methodologies. The instantaneous SNR is defined by SNR= 
$ 3 (w)/* v (u;). For PS, the optimal subtraction factor 6 is given 
by (58) and for IPS, G(w) is given by (45) with $xM and $ v (u>) 
there replaced by, respectively, $ z (u/) and & v [u>). 

H{u>) Bias Variance 

E[*,(w)]/» tf (w) Var(* a (w))/ 7 * 2 (u;) 

SPS 1-6 . 6 2 

MS -2(v/l + SNR-l) (v/l+SNR- l) 2 

;pc ^SNR ( SNR 2 \ 2 f, « i+SNR V 

IPS "t+SNR 7 ISNR 7 ^; V + 7 SNI^J 

wf _SNR 4 ( SNR \ 2 

A/L ^(v'SNR+l- v'SNR) 2 ^ (l + ^1 + 
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p speech-in-noise model order 

p running average update factor for $„(w) 

3. For each frame of input data do: 

(a) Speech detection (step 110) 

The variable Speech is set to true if the VAD output equals st = 21 or st = 22. 
Speech is set to false if st = 20. If the VAD output equals st = 0 then the 
algorithm is reinitialized. 

(b) Spectral estimation 

If Speech estimate $ z (u;): 

i. Estimate the coefficients (the polynomial coefficients {A,. . Jp} and the 
variance a*) of the all-pole model (23) using the autocorrelation method 
applied to zero mean adjusted input data {x{k)} (step 120). 

ii. Calculate * x {u) according to (25) (step 130). 

else estimate 4> v (w) ( ste P 14 °) 
i. Update the background noise spectral model *„(w) using (4), where 
is the Periodogram based on zero mean adjusted and Hanning/Hamming 
windowed input data x. Since windowed data is used here, while * r (w) 
is based on unwindowed data, *„(«) has to be properly normalized. A 
suitable initial value of # w (w) is given by the average (over the frequency 
bins) of the Periodogram of the first frame scaled by, for example, a factor 
0.25, meaning that, initially, a apriori white noise assumption is imposed 
on the background noise, 
(c) Spectral subtraction (step 150) 

i. Calculate the frequency weighting function H (w) according to Table 1. 

ii. Possible postfiltering, muting and noise floor adjustment. 

iii. Calculate the output using (3) and zero-mean adjusted data {x(k)). The 
data {x(Jfc)} may be windowed or not, depending on the actual frame 
overlap (rectangular window is used for non-overlapping frames, while a 
Hanning window is used with a 50% overlap). 
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From the above description it is clear that the present invention results in a sig- 
nificant noise reduction without sacrificing audible quality. This improvement may be 
explained by the separate power spectrum estimation methods used for speech and non- 
speech frames. These methods take advantage of the different characters of speech and 
non-speech (background noise) signals to minimize the variance of the respective power 
spectrum estimates 

• For non-speech frames $ v (u/) is estimated by a non-parametric power spectrum 
estimation method, for example an FFT based periodogram estimation, which uses 
all the N samples of each frame. By retaining all the N degrees of freedom of the 
non-speech frame a larger variety of background noises may be modeled. Since the 
background noise is assumed to be stationary over several frames, a reduction of the 
variance of $ v (u>) may be obtained by averaging the power spectrum estimate over 
several non-speech frames. 

• For speech frames $ x {u) is estimated by a parametric power spectrum estimation 
method based on a parametric model of speech. In this case the special character 
of speech is used to reduce the number of degrees of freedom (to the number of 
parameters in the parametric model) of the speech frame. A model based on fewer 
parameters reduces the variance of the power spectrum estimate. This approach is 
preferred for speech frames, since speech is assumed to be stationary only over a 
frame. 

It will be understood by those skilled in the art that various modifications and changes 
may be made to the present invention without departure from the spirit and scope thereof, 
which is defined by the appended claims. 
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ANALYSIS OF #msM 

Paralleling the calculations for Hms{u) gives 

2 



*.(")= 1 



1- 



N 



(27) 



where in the second equality, also the Taylor series expansion Vl+x ~ 1 + x/2 is used. 
From (27) it follows that the expected value of $,(u/) is non-zero, given by 



£[*,(«)] ^ 2* v («") 1 - 



(28) 



Further, 



Var(*,(u>)) ^ 



(29) 



1-. 



Var(*,(w)) + Var(* w (u/)) 



Combining (29) and (15) 

Var{$ 3 (")) - [l - 



1 + 



♦•M 1 -. *2 



7*>) 



(30) 
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APPENDIX B 

ANALYSIS OF Hwf(") 

In this Appendix, the PSD error is derived for speech enhancement based on Wiener 
filtering, [2], In this case, H{u) is given by 

Hwf{") = - , *' (< "? , , = Hpsi") (31) 

Here, <J> 3 (u) is an estimate of $ 3 {u>) and the second equality follows from $,(u/) = $ x (u;) — 
4 v (u). Noting that 



a straightforward calculation gives 



K"* vH + 2 {^M Ax(w) " Au(u;) }) 



From (33), it follows that. 



and 

2 



(32) 



(33) 



E[i.M)*-(l-%^* v {u) . (34) 



Var(».(«))=f4(l-|^jJ 7 *2(") (35) 
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APPENDIX C 



ANALYSIS OF H ML {w) 

Characterizing the speech by a deterministic wave-form of unknown amplitude and 
phase, a maximum likelihood (ML) spectral subtraction method is defined by 



H ML (u>) = - \ 1 + 



= i(i + tfpsM) 

Inserting (11) into (36) a straightforward calculation gives 



(36) 



)A x ( w )Vl 



* 2 1 + 



(37) 



1 1 / ^v(^) A / \ A / \ 



where in the first equality the Taylor series expansion (1 + x) 1 ^ 1 - x and in the 
second y/T+x ~ 1 + x/2 are used. Now, it is straightforward to calculate the PSD 
error. Inserting (37) into (9)-(l0) gives, neglecting higher than first order deviations in 
the expansion of Hml{u) 



<*>,(")- -(1 + 



4> x (w) - $ 3 {uj) 



+ -4 1 + 



A z (u/) - A v (lj) 



) 



(38) 



From (38), it follows that 



4> x (u;) - $ 3 {u) 



(39) 



where in the second equality (2) is used. Further, 



Va,(*,M).-li + .-M- w 



7*:(w) 



(40) 
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APPENDIX D 

DERIVATION OF H iFS {v) 

When 4> x (u;) and 4> v (u;) are exactly known, the squared PSD error is minimized by 
Hpsiv), that is Hps{v) with $ x {u) and 4 v (a;) replaced by * x (u>) and * v (w), respectively. 
This fact follows directly from (9) and (10), viz. $ s (u) = [J/ 2 (u/)* x (ui) - $ 5 (^)] 2 = 0 : 
where (2) is used in the last equality. Note that in this case H(u) is a deterministic quan- 
tity, while H{u) is a stochastic quantity. Taking the uncertainty of the PSD estimates into 
account, this fact, in general, no longer holds true and in this Section a data-independent 
weighting function is derived in order to improve the performance of Hp${u)). Towards 
this end, a variance expression of the form 

Var(* 5 H)sf 7 »;H (41) 

is considered (£ = 1 for PS and £ = (1 - \/l + SNR) 2 for MS and 7 = 7* + 7«)- The 
variable 7 depends only on the PSD estimation method used and cannot be affected by 
the choice of transfer function H{u). The first factor however, depends on the choice 
of J/(u/). In this section, a data independent weighting function G(u) is sought, such that 
H(u) = \jG(u) Hps{uj) minimizes the expectation of the squared PSD error, that is 

G{u) = argminE[*3(u;)] 2 

(42) 

= G(«) H 2 PS {u>) i> x (u) - 9 a (u>) 

In (42), G(u) is a generic weigthing function.. Before we continue, note that if the weight- 
ing function G{u) is allowed to be data dependent a general class of spectral subtraction 
techniques results, which includes as special cases many of the commonly used methods, 
for example, Magnitude Subtraction using G(u/) = H 2 ms {uj)/ H 2 ps {uj). This observation 
is, however, of little interest since the optimization of (42) with a data dependent G(w) 
heavily depends on the form of G(u/). Thus the methods which use a data-dependent 
weighting function should be analyzed one-by-one, since no general results can be derived 
in such a case. 

In order to minimize (42), a straightforward calculation gives 
*,(w) = (C(u)- !)*,(") 
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<43) 

Taking expectation of the squared PSD error and using (41) gives 

E[*,(u)) 7 * - l) a *2(") + <? 2 (u/) 7 *» (44) 

Equation (44) is quadratic in G(w) and can be analytically minimized. The result reads, 

G(w) ~ •;.(«) +7*50*) (45) 

1 

~ i i , r 

where in the second equality (2) is used. Not surprisingly, 6{u>) depends on the (unknown) 
PSDs and the variable 7. As noted above, one cannot directly replace the unknown 
PSDs in (45) with the corresponding estimates and claim that the resulting modified PS 
method is optimal, that is minimizes (42). However, it can be expected that, taking the 
uncertainty of * x (w) and *.(«) into account in the design procedure, the modified PS 
method will perform "better" than standard PS. Due to the above consideration, this 
modified PS method is denoted by Improved Power Subtraction (IPS). Before the IPS 
method is analyzed in APPENDIX E r the following remarks are in order. 

For high instantaneous SNR (for u such that *,(w)/*vM » 1) it follows from (45) 
that. G(u») ~ 1 and, since the normalized error variance Var(* # (u/))/*2(w), see ( 41 ) is 
small in this case, it can be concluded that the performance of IPS is (very) close to the 
performance of the standard PS. On the other hand, for low instantaneous SNR (for w 
such that 7 *2M » <5(w) * *»/(7* 2 »). fading to, cf. (43) 

and 

w-w-M (47) 

However, in the low SNR it cannot be concluded that (46)-(47) are even approximately 
valid when G(u) in (45) is replaced by that is replacing <*>*(") and *„(<") in (45) 
• with their estimated values and respectively. 
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APPENDIX E 

ANALYSIS OF H IPS (w) 

In this APPENDIX, the IPS method is analyzed. In view of (45), let G(lj) be defined 
by (45), with $„(u>) and ^(w) there replaced by the corresponding estimated quantities. 
It may be shown that. 

* 4 (u>) s (G(u.O - l)* a (w) 

(llR Ax(w) " AvH ) (48) 

which can be compared with (43). Explicitly, 

E[*.{u>)) ^ (<5(w) ~ 1)*» 

and 

Var(*,(u/)) as G 2 (w) 

For high SNR, such that * a (u;)/*„(u) » 1, some insight can be gained into (49)-(50). In 
this case, one can show that 

■£(*,(«)] =5-0 (51) 

and 

Var($ 4 (u;)) = + ^rf^) 7^M (52) 

The neglected terms in (51) and (52) are of order 0((* v (u;)/* a (w)) 2 ). Thus, as al- 
ready claimed, the performance of IPS is similar to the performance of the PS at high 
SNR. On the other hand, for low SNR (for u such that *5(w)/(7*2(w)) < - 
*2(«)/(7*2M). and 

£[*,(«)] = -*,(«) (53) 



(49) 



(50) 
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and 

v«(*»)*»^) (54) 

Comparing (53)-(54) with the corresponding PS results (13) and (16). it is seen that 
for low instantaneous SNR the IPS method significantly decrease the variance of •.(«) 
compared to the standard PS method by forcing *.(«) in (9) towards zero. Explicitly, 
the ratio between the IPS and PS variances are of order 0(*<»/* 4 »). One may also 
compare (53)-(54) with the approximative expression (47), noting that the ratio between 
them equals 9. 



2ti 
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APPENDIX F 

PS WITH OPTIMAL SUBTRACTION FACTOR 6 

An often considered modification of the Power Subtraction method is to consider 



\ *x(w) 



where 6(uj) is a possibly frequency dependent function. In particular, with 6(u) = 6 for 
some constant S > 1, the method is often referred as Power Subtraction with oversub- 
traction. This modification significantly decreases the noise level and reduces the tonal 
artifacts. In addition, it significantly distorts the speech, which makes this modification 
useless for high quality speech enhancement. This fact is easily seen from (55) when 6 > 1. 
Thus, for moderate and low speech to noise ratios (in the u>-domain) the expression under 
the root-sign is very often negative and the rectifying device will therefore set it to zero 
(half- wave rectification), which implies that only frequency bands where the SNR is high 
will appear in the output signal s(k) in (3). Due to the non-linear rectifying device the 
present analysis technique is not directly applicable in this case, and since 6 > 1 leads to 
an output with poor audible quality this modification is not further studied. 

However, an interesting case is when 6{u) < 1, which is seen from the following 
heuristical discussion. As stated previously, when $ x (u) and $ v (u>) are exactly known, 
(55) with 6{u) = 1 is optimal in the sence of minimizing the squared PSD error. On the 
other hand, when $ x (u/) and 4> v (u;) are completely unknown, that is no estimates of them 
are available, the best one can do is to estimate the speech by the noisy measurement 
itself, that is s(k) = x(fc), corresponding to the use of (55) with 6 = 0. Due the above 
two extremes, one can expect that when the unknown $ x (u) and $ v (u) are replaced by, 
respectively, $ x (v) and 4> v (w), the error E[$ 3 (lj)) 2 is minimized for some 6(u) in the 
interval 0 < 8(u) < 1. 

In addition, in an empirical quantity, the averaged spectral distortion improvement, 
similar to the PSD error was experimentally studied with respect to the subtraction factor 
for MS. Based on several experiments, it was concluded that the optimal subtraction factor 
preferably should be in the interval that span from 0.5 to 0.9. 
Explicitly, calculating the PSD error in this case gives 

* (i - «H)* w (o/) + «h ^M^H - A v (u) 
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(56) 



Taking the expectation of the squared PSD error gives 

El*.(u,)] a = (1 - 6{»)?**M + * 8 ^) (5?) 
where (41) is used. Equation (57) is quadratic in «(«) and can be analytically minimized. 
Denoting the optimal value by 5, the result reads 

J". — < 1 < 58) 
1+7 

Note that since 7 in (58) is approximately frequency independent (at least for N » 1) 
also 6 is independent of the frequency. In particular, 6 is independent of ♦.<„) and •„(«>. 
which implies that the variance and the bias of •.(*) directly follows from (57). 

The value of I may be considerably smaller than one in some (realistic) cases. For 
example, once again considering , v = 1/r and 7 , = 1. Then 6 is given by 

- 1 1 
6 = -■ 



2 1 + 1/2t 

which, dearly, for an , is Smaller than 0.5. .» this case, the fact that i « 1 indicates 
that the uncertainty in the PSD estimators (and, in particular, the uncertainty in »,M) 
have a .arge impact on the quality (in terms of PSD error) of the output. Especially, the 
use of I « 1 implies that the speech to noise ratio improvement, from input to output 

signals : is small. 

An arismg question is that if there, similarly to the weighting function for the IPS 
method in APPENDIX D, exists a data independent weighting function G(w) In AP- 
PENDIX G, such a method is derived (and denoted HPS). 
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APPENDIX G 

DERIVATION OF H S!PS (u>) 

In this appendix, we seek a data independent weighting factor C(u>) such that H{u) = 
y/^)H S psH for some constant 6 (0 < 6 < 1) minimizes the expectation of the squared 
PSD error, cf (42). A straightforward calculation gives 



*» = (C(u>) - l)*,(u») + G(u;)(l - 6)* v (u) 

The expectation of the squared PSD error is given by 

£[*,M] 2 = (G(u,) - 1) 2 $J( W )+G 2 ( U )(1 - 6) 2 *l(u) 

2(G(w) - l)*,(w)C(w)(l - *)* w (w)+C 2 (u;)* 2 7*2(w) 



(59) 



(60) 



The right hand side of (60) is quadratic in G(u) and can be analytically minimized. The 
result G(u>) is given by 

G(u,) = 

tfM + *.(h/)*«,(h*)(l-f) 



*2(w) + M J (w)* v (u/)(l-«) + (l-i)2*2( u/ ) + *2 7 *2( w ) 

1 

(61) 



^ I ♦,(*)-*«(«-) J 

where 0 in the second equality is given by 

(1 - S) 2 + 6 2 <y + (1 - f )».(u;)/» v (w) 

1 + (1 - «)*«(«)/*.(«) . (62) 

For (5 = 1, (61)-(62) above reduce to the IPS method, (45), and for 6 = 0 we end up 
with the standard PS. Replacing *,(«) and in (61)-(62) with their corresponding 
estimated quantities * x ( w ) - and *„( w ), respectively, give rise to a method, which 
in view of the IPS method, is denoted SIPS. The analysis of the f IPS method is similar to 
the analysis of the IPS method, but requires a lot of efforts and tedious straightforward 
calculations, and is therefore omitted. 
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CLAIMS 

1. A spectral subtraction noise suppression method in a frame based digital commu- 
nication system, each frame including a predermined number N of audio samples, thereby 
giving each frame N degrees of freedom, wherein a spectral subtraction function H(uj) is 
based on an estimate $ v (u) of the power spectral density of background noise of non- 
speech frames and an estimate $ x (u) of the power spectral density of speech frames, 
characterized by: 

approximating each speech frame by a parametric model that reduces the number of 
degrees of freedom to less than N; and 

estimating said estimate $ x {u>) of the power spectral density of each speech frame by 
a parametric power spectrum estimation method based on the approximative parametric 
model 

estimating said estimate $ v {uj) of the power spectral density of each non-speech frame 
by a non-parametric power spectrum estimation method. 

2. The method of claim 1, characterized by said approximative parametric model 
being an autoregressive (AR) model. 

3. The method of claim 2, characterized by said autoregressive (AR) model being 
approximately of order 

4. The method of claim 3, characterized by said autoregressive (AR) model being 
approximately of order 10. 

5. The method of claim 3, characterized by a spectral subtraction function H{uj) in 
accordance with the formula: 



//(u/) = 



1 



z(w)/ 



where G(u>) is a weighting function and 8(lj) is a subtraction factor. 1 

6. The method of claim 5, characterized by G{u) — L 

7. The method of claim 5 or 6, characterized by 6{u) being a constant < 1. 

8. The method of claim 3, characterized by a spectral subtraction function H(u>) in 
accordance with the formula: 

Hlu) = 1 - * 
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9. The method of claim 3, characterized by a spectral subtraction function H{u) 
accordance with the formula: 



■ 10. The method of claim 3, characterized by a spectral subtraction function H(u>) 
accordance with the formula: 



^ 
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